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AiResearch  Manufacturing  Company  of  Arizona  has  conducted 
a research  study  for  the  Office  of  Naval  Research  under  Contract 
N00014-74-C-0 317 , Task  NR  061-221,  to  develop  rapid  efficient 
techniques  for  Navier-Stokes  solutions.  This  document  reports 
the  last  year's  activity.  Previous  activities  have  been 
reported  and  are  referenced  herein. 
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SECTION  I 
INTRODUCTION 


Traditional  techniques  of  solving  the  Navier-Stokes  equa- 
tions for  compressible  separated  flow  have  utilized  time- 
dependent  methods.  These  methods  have,  with  extensive  develop- 
ment, steadily  improved  in  both  accuracy  and  computational  times. 
Recently,  however,  Dodge^  and  Dodge  and  Lieber^  have  introduced 
a new  steady-state  solution  method  for  incompressible  flow  that 
splits  the  pressure  and  shear  forces  into  two  separate  equations. 
One  equation  is  elliptic,  and  the  other  nearly  parabolic.  Con- 
tained herein  is  a description  of  the  first  attempt  to  extend 
this  method  to  compressible  flow.  Application  is  made  to  the 
interaction  of  an  oblique  shock  wave  with  a laminar  boundary 
layer.  The  development  was  guided  by  the  realization  that  the 
ultimate  value  of  such  a method  lies  not  in  the  solution  of  two- 
dimensional  laminar  shock  boundary  layer  interactions,  but  rather 
in  applications  to  more  complete  two-dimensional  and  three- 
dimensional  geometries.  Although  more  development  to  improve 
accuracy  and  further  reduce  computational  times  is  suggested, 
the  results  to  date  demonstrate  the  general  viability  of  this 
numerical  method.  The  following  sections  describe  in  detail  the 
basic  equation  development,  numerical  techniques,  results  of  the 
calculations,  and  suggestions  for  accuracy  improvements. 

SECTION  II 

EQUATION  DEVELOPMENT 


The  basis  of  Dodge's  equation-splitting  solution  method 
lies  in  the  separation  of  equations  describing  pressure  and 
viscous^ef fects . To  accomplish  this  separation,  the  velocity 
vector  W is  expressed  in  the  following  manner: 

W = V0  + U (1) 

where  <t>  = velocity  potential 

U = velocity  vector  reflecting  viscous  effects 
(viscous  loss  velocity  vector) 

Equations  are  then  required  to  solve  individually  for  4>  and 
U. 

As  an  initial  step,  consider  the  vector  form  of  the  momen- 
tum equation  for  steady  laminar  flow: 


p (W*V) W 


VP  + V(XV*W) 


+ V • (pdefW) 


(2) 


/ 


A 


1 
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where  p = density 

P = static  pressure 

u = molecular  viscosity 

X = second  coefficient  of 
viscosity 

defW  = VW  + (VW)* 

(VW)*  = transpose  (VW) 

Equation  (1)  may  then  be  substituted  for  W in  Equation  (2)  as 
follows : 

p V<J)  • V ( VcJ> ) + pU*V(V<J>)  + p(W*V)U 
= - VP  + V(AV*W)  + V • ( pdefW) 


To  treat  the  static  pressure  gradient  in  Equation  (3),  consider 
the  density  to  be  described  by  two  terms  — one,  p*,  related 
to  static  pressure  and  the  other,  py,  related  to  viscous 
effects: 

P = P*  + Pv  (4) 


The  density  component  p*  is  defined  by 


where 


P* 


p;a 


i 

Y-l 


Pj  = reference  condition  of  stagnation  density 
Tj  = reference  condition  of  stagnation  temperature 
Cp  = specific  heat  at  constant  pressure 
y = ratio  of  specific  heats 


(5) 
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Further,  one  can  show  that  all  possible  static  pressure  gradi 
ents  may  be  represented  by 

- 

1 

- VP  = p*  V$*V(V$) 

(6) 

1 

1 

where  p*  is  selected  so  that  the  right-hand  side  of  Equation 
is  irrotational . Substituting  Equation  (4)  into  Equation  (3) 
combining  terms,  and  using  Equation  (6)  yields 

(6) 

9 

1 

(p  - p*)  V<J>*V  (V$>)  + pU*V(V<(>)  + p(W*V)U 

(7) 

1 

= V (XV*ft)  + V* (pdefW) 

f 

r 

I 

Equations  (6)  and  (7)  co^ld  then  provide  the  basis  for 
separate  equations  for  <J>  and  U,  respectively.  However,  when 
Equation  (6)  is  combined  with  the  continuity  equation  and  the 
equation  of  state,  it  yields  a relation  that  becomes  hyperbolic 
for  supersonic  conditions  based  upon  potential  rather  than 
total  velocity.  To  reconcile  this  difficulty,  a different 
approach  was  taken  in  the  development  of  the  potential  equation. 

The  pressure  gradient  may  be  expressed  in  terms  of  the 
equation  of  state,  as 

** 

VP  = a2Vp  + a2p(I^)VS 

(8) 

r 

where 

a = local  speed  of  sound  at  static  conditions 

u 

R = gas  constant 

S = entropy 

| 

f 

Replacing  the  pressure  gradient  term  in  the  momentum  equation 
Equation  (2),  with  Equation  (8),  and  dotting  the  resulting 
expression  with  ft,  yields  Equation  (9) . 

9 

t 

ft*  [p(W*V)W]  = - a 2 (ft*  Vp ) - a2p  (3-~-1-)  W«VS 

Y K 

+ W-V(XV*W)  + ft* [V* (pdefW) ] 

(9) 

f 

The  vector  form  of  the  compressible  continuity  equation,  given 
by 

F 

m 

o 

II 

+9 

Q. 

• 

> 

(10) 

I 

1 

3 

4 

||SE 

may  be  expanded  and  re-expressed  as 


w* Vd  = - pv  • w 


(11) 


Substituting  Equation  (11)  into  the  first  term  of  the  right-hand 
side  of  Equation  (9),  rearranging,  and  dividing  by  p,  yields 
Equation  (12). 


W- [ (W*7) W]  - a (V*W) 


= - a2(I 


yR 


)W*  VS 


+ w • [ —V ( X V • W ) ] + W*  [-V*  (ydefW) ] 

P P 


(12) 


However,  Equation  (12)  may  be  recast  by  replacing  U in  the  left- 
hand  terms  of  Equation  (12)  with  the  velocity  components  of 
Equation  (1).  Then,  placing  all  terms  involving  potential  on 
the  left-hand  side  of  the  resulting  equation,  and  expanding  the 
last  two  terms,  yields  the  following: 


W*[(W*V)V<t>]  - a2v2(j>  = a 2v  • u - W*[(W*V)U] 
-a  2 (^y"--1-)  W- VS  + w.  [iv(XV*W)]  (13) 

+ W*{ji[V(V*W)  + V2W]  + (defW)^H-} 


Equation  (13)  thus  describes  pressure  effects  in  terms  of  the 
velocity  potential,  <p.  The  equation  will  be  elliptic  or  hyper- 
bolic, depending  upon  the  local  Mach  number,  defined  in  terms  of 
the  velocity,  W. 

The  equation  for  the  viscous  loss  velocity  vector  is 
developed  from  the  momentum  equation,  as  presented  in  Equation 
(7) . The  terms  on  the  right-hand  side  of  Equation  (7)  may  be 
expanded  in  the  same  manner  as  for  the  incompressible  case  pre- 
sented by  Dodge1.  Then,  rearranging,  and  dividing  the  resulting 
equation  by  p,  yields 

U*  V (V<p)  + (W-V)U  = (£p  - 1)  V<t>*  V (V<f>) 

+ i-  { ( v2  4> ) VX  + (V-U)VX  + ( 2p  + X)  V ( V 2 <J> ) 

+ (p  + X)V(V*U)  + uV2U  + 2 [ V ( V<{> ) ] * Vp 


(defU) Vp} 


+ 


(14) 


When  expanded.  Equation  (14)  will  produce  scalar  equations 
for  the  viscous  loss  velocity  components,  U^.  These  equations 
may  be  parabolized  by  treating  the  streamwise  diffusion  term  as 
a known  quantity  from  the  previous  iteration,  as  discussed  by 
Dodgel  and  Dodge  and  Lieber^. 

Because  the  flow  is  compressible,  an  energy  equation  is 
also  required.  In  terms  of  static  enthalpy,  the  laminar  steady- 
flow  energy  equation  may  be  given  as 

p (W*  V)  (h  + iw2)  = V* (XWV*W) 

- - (15) 

+ V*(W*pdefW)  + V*(k7T) 

where 

h = enthalpy  at  static  conditions 
k = conductivity  of  the  fluid 
T = static  temperature 

Expressing  enthalpy  in  terms  of  static  temperature,  i.e.. 


h = C T (16) 

P 


and  expanding  the  left-hand  side  of  Equation  (15),  yields 
Equation  (17) 


p (W*  7)  (CpT)  + i p (W-V)W2  = 
+ V* (W*  pdefW)  + V*  (kVT) 


V* (XW7*W) 

(17) 


Assuming  a constant  Prandtl  number,  Pr,  given  by 


Pr 


(18) 


then  conductivity  may  be  expressed  as 
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In  addition,  specific  heat  is  assumed  to  be  constant.  Thus, 
substituting  Equation  (19)  into  Equation  (17),  rearranging,  and 
dividing  by  pCp,  yields 

(W-V)T  - -1-  y.  [ (^R)VTJ 

= [“  (W*V)W2  (20) 

+ V*  ( AWV • W)  + v • (W*  ydefW) ] 

Expanding  the  second  term  on  the  left-hand  side  of  Equation  (20) 
results  in  the  following: 


(W-V)T  - — ^ [VyVT  + uV2T] 

= ■—  [-  j (W-V)W2  + V • ( AWV  • W) 

P P 

+ V*  (W- ydefW) ] (21) 


Expanding  the  last  two  terms  of  Equation  (21)  yields 

(W- V ) T - (Vu-VT)  - ^ (uV2T) 

= ~~  {-  j (W-V)W2  + A (V-W)  2 
P 

+ U [ V • (W  • defW)  - W- ( V-defW) ] 

+ W* [V ( X V • W)  + V • ( udefW) ] } (22) 


Substituting  the  momentum  equation  in  the  form  of  Equation  (7) 
for  the  last  two  terms  in  Equation  (22)  gives  the  final  form  of 
*-he  laminar,  steady-flow  energy  equation  in  terms  of  static 
temperature : 

(W-V)T  - ■—  (Vu-VT)  - (uV2T) 


1 

PC. 


1 2 


(W-V)W  + A (V-W) 


+ u [V-(W-defW)  - W-(V-defW)) 


+ pW-  [U-V(V«j»)  + (W-V)U 
+ (1  - jp)  V<D - V (V<|>)  ] } (23) 


6 


4 


I 


I 

I 


iT 


y 


637-7° 


As  with  the  equation  for  U,  Equation  (23)  may  be  parabolized  by 
treating  the  diffusion  term  as  a known  quantity  from  the  previous 
iteration. 

Equations  (13),  (14),  and  (23)  comprise  the  system  of  equa- 

tions that  must  be  solved  to  yield  <J> , U,  and  T.  In  conjunction 
with  these  equations,  a set  of  boundary  conditions  is  required. 
Figures  1 and  2 illustrate  the  boundary  conditions  on  the  solu- 
tion space  for  the  elliptic/hyperbolic  and  parabolic  equations. 

In  Figure  1,  the  potential  velocities  Vj , VM,  and  VE  are 
determined,  based  upon  the  oblique  shock  equations.  The  value 
of  is  determined  by  integrating  the  values  of  9<()/9x^  to  the 
downstream  boundary. 

The  upstream  values  of  U ancj  T in  Figure  2 are  dependent 
upon  the  corresponding  value  of  W,  which  is  initialized  using  an 
arbitrary  boundary  l^iyer  profile  similar  to  a Pohlhausen  profile. 
An  isothermal  plate  is  assumed,  with  plate  temperature  equal  to 
the  reference  stagnation  temperature. 


7 
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SECTION  III 
NUMERICAL  METHODS 

Two  major  numerical  solutions  are  required  to  implement  the 
above  algorithm:  a solution  to  a mixed  elliptic/hyperbolic 

potential  equation,  and  a solution  to  a nearly  parabolic  viscous- 
shear  equation.  Murman  and  Cole^,  in  a pioneering  paper, 
developed  a technique  fo'"  modifying  a difference  star  with  Mach 
number,  and  then  using  conventional  relaxation  to  solve  a 
similar  equation.  Modifications  to  the  Murman  and  Cole  method 
have  been  necessary  to  improve  its  accuracy,  and  allow  its  use 
in  totally  supersonic  flow.  Dodge'*  *5  describes  such  a method 
and  its  application  to  supersonic  inlet,  subsonic  exit  compressor 
blades.  Murman® , in  a subsequent  paper,  reported  further  improve- 
ments in  the  basic  methods  to  account  for  differencing  across 
discontinuities.  From  these  works,  a list  of  ideal  characteris- 
tics can  be  assembled  for  the  supersonic  portion  of  any  such 
algorithm. 

o Region  of  influence  of  difference  star  matches 
precisely  the  region  of  influence  of  the 
differential  equation. 

o Across  a potential  slope  discontinuity,  the 

method  must  preserve  the  proper  jump  conditions. 

o The  method  remains  integrally  conservative  in  the 
non-discontinuous  portion  of  the  flow  field. 

Unfortunately,  such  a numerical  method  does  not  now  exist. 

To  assist  in  the  evaluation  of  the  importance  of  each  of 
the  above  criteria,  the  one-dimensional  potential  Equation  (24) 
was  studied. 


2,d‘ 
a ) — 


dx 


= 0 


(24) 


where 


a 

V 


sound  speed 

d <t> 

d£ 


= velocity  potential 


9 


s~  . 
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For  the  solution  to  be  other  than  trivial,  the  flow  must  enter 
supersonically,  undergo  a jump,  and  exit  subsonically . Boundary 
conditions  then  consist  of  an  inlet  Mach  number  and  an  exit 
potential.  By  integrating  Equation  (24)  and  assumina  an 
ideal  gas,  a jump  condition  can  be  derived  that  relates  down- 
stream velocity  to  inlet  velocity. 


■ - v v > 

(a'2  - Y-  g.  1 V22)  (25) 


where 

a'  = stagnation  speed  of  sound 
•y  = ratio  of  specific  heats 
= velocity  ahead  of  the  jump 
V?  = velocity  behind  the  jump 

The  analytical  solution  ahead  of  and  behind  the  jump  is  obviously 
that  of  constant  velocity.  Jump  location  is  determined  by 
matching  to  the  exit  potential  boundary  condition  using 
Equation  (26) . 


XJ  = 


4>  - V„L 

^EXIT  2 

V - V 
1 2 


(26) 


where 

xj  = position  of  the  jump 
L = total  length  of  the  solution 

region  (X  = 0 is  the  inlet) 

<J>Exit  = potential 

Consider  what  happens  when  the  original  Murman  and  Cole  method 
is  applied.  Upstream  of  the  assumed  jump  location,  the  differ- 
ence star  is  cast  in  a backwards  direction. 

(<Jk  - 24k_^  + ) = 0 (27) 


where 


4k  = potential  at  the  ith  node 


10 
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Note  that  the  value  of  - ar  is  immaterial  since  it  is  a 
non-zero  constant  at  any  particular  node,  and  thus  divides  out. 
Downstream  of  the  assumed  jump  location,  the  equation  is 
elliptic. 


*i+l  ' 2*i  * *i-l  = ° 


The  solution  to  both  Equations  (27)  and  (28)  is  a linear  varia- 
tion in  $ . An  initial  guess  Tor  potential  would  be  made. 

Equation  (27)  is  solved  starting  at  node  3,  and  simply  broadcasts 
forward  the  initial  velocity  gradient  until  a node  is  reached, 
which  was  initially  assumed  subsonic.  At  this  point,  the  sub- 
sonic solution  proceeds  as  a boundary  value  problem  with  <p  at 
the  change-over  (supersonic  to  subsonic)  point  established  by 
Equation  (27).  There  is,  however,  no  communication  between 
upstream  and  downstream.  Thus,  the  location  of  the  jump  never 
changes,  and  the  value  of  the  velocity  at  the  jump  is  only  con- 
trolled by  the  initial  guess  of  jump  location  and  potential 
variation.  No  aspect  of  the  physical  problem  appears  in  the 
numerical  solution  at  the  jump.  To  overcome  this,  special  condi- 
tions at  the  jump  must  be  invoked,  as  pointed  out  by  Murman^. 

It  is,  however,  immaterial  what  form  is  used  away  from  the  jump. 

A conservative  form  results  in  a solution  that  is  no  different 
from  a nonconservative  form. 

When  actually  applying  methods  to  introduce  upstream- 
downstream  communication,  two  classes  result.  The  first  class, 
into  which  Wurman ' s method  falls,  allows  the  downstream  flow  to 
be  adjusted  based  on  upstream  conditions.  This  class  allows 
only  for  jump  propagation  in  the  forward  direction.  The  second 
class,  permits  the  jump  to  retreat  by  allowing  the  upstream  flow 
to  be  affected  by  downstream  conditions.  Whether  the  jump  needs 
to  retreat  or  go  forward  is  only  dependent  on  the  initial  guess, 
and  thus,  for  complete  generality,  both  classes  must  be  combined 
into  a single  method.  For  one-dimensional  flow,  such  techniques 
are  readily  generated,  the  results  of  which  are  shown  in 
Figure  3 for  a forward  propagating  jump  and  Figure  4 for  a 
retreating  jump.  The  jump-node  region  is  treated  by  the  fol- 
lowing algorithm 


D = r<p . , + 4> . , - (1+r)  <j>  • = 0 
l+l  i-l  i 


applied  to  the  node  ahead 


of  the  jump  if 


D >0, 
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Figure  3. 


Forward  propagating  one-dimensional  jump. 
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applied  to  the  node  behind  the  jump  if  D <"0 


aDOWN 


a = V (a  ' 2 


v2) 


where  DOWN  and  UP  subscripts  refer 

to  downstream  and  upstream  of  the 

jump. 

It  is  necessary  to  evaluate  the  velocity  away  from  the  jump, 
since  values  at  the  jump,  particularly  for  intermediate  ite- 
rations, are  in  significant  error. 

When  the  flow  field  is  multidimensional  and  jumps  are  other 
than  normal,  the  proper  procedure  is  complicated  by  the  lack  of 
accepted  numerical  solutions.  Murman's  method  assumes  an 
oblique  jump  passing  through  the  grid  system  at  an  arbitrary 
angle.  However,  this  derivation  looses  sight  of  the  fact  that 
the  differential  equation  is  left  behind  when  solutions  are  made 
to  the  finite-difference  grid.  Discontinuities  on  a fixed 
orthogonal  finite-difference  system  can  only  have  one  of  three 
forms,  illustrated  by  Figure  5.  A physically  oblique  wave  must 
be  made  up  of  discrete  components  of  the  three  jumps.  The  normal 
jump  causes  little  or  no  problems,  but  the  other  two  result  in 
differencing  across  a discontinuity.  To  be  totally  correct,  as 
Murraan  points  out,  any  difference  scheme  applied  to  any  combina- 
tion of  the  above  discontinuities  must  be  integrally  (summation 
in  difference  system)  conservative  of  the  desired  jump  condition. 
Even  if  such  a system  existed,  however,  practical  problems,  such 
as  detecting  a jump  and  determining  proper  upstream  and  down- 
stream conditions  away  from  the  perturbed  region  of  the  jump, 
intervene  to  prevent  reliable  treatment  of  jump  conditions. 

The  method  by  Dodge**  avoids  some  of  the  difficulty  by  only 
differencing  along  characteristic  lines,  and  not  differencing 
across  oblique  jumps,  which  may  explain  some  of  its  good  results. 
However,  it  requires  a nonstationary  grid  system  that  is  diffi- 
cult to  implement  in  two  dimensions,  and  is  nearly  impossible  in 
three  dimensions.  Thus,  after  comparing  other  methods,  such  as 
the  two  by  Murman,  to  a range  of  six  supersonic  and  supercritical 
cases,  it  was  decided  to  mix  a nonorthogonal  and  orthogonal 
difference  system,  staying  with  a fixed  orthogonal  grid  system. 
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The  deqree  of  mixing  was  determined  as  a function  of  local 
Mach  number,  based  upon  a criterion  developed  from  a stability 
analysis  of  the  basic  relaxation  difference  equation.  Figure  6 
illustrates  the  difference  stars  utilized.  In  simplified  form, 
the  mixing  of  difference  operators  is  accomplished  in  the  following 
following  manner: 


r s24'  + r ‘)2^-  - re  9-2-^  + C — -]  (f  ) 
C1  ^ 2 9y2  ' [ 1 9X2  2 9y2  H 


+ [Cl  4 + c2 

1 9x^  2 3y 


- fH> 


where  f„  represents  the  fraction  of  hyperbolic  difference  opera- 
tor utilized,  and  (1-fj.j)  represents  the  fraction  of  single-spaced 
implicit  difference  operator.  The  local  value  of  f^  was  deter- 
mined for  supersonic  flow  to  be 


f H = -C2/C1  for  (~C2/C1)  < 1.0 

or 

fH  = 1.0  for  {-C2/C\)  i 1.0 

based  upon  a stability  analysis. 


If  the  variable  supersonic  difference  star  is  examined  in 
relation  to  the  location  of  local  characteristic  lines,  it  may 
be  seen  that,  as  Mach  number  increases,  the  difference  star 
becomes  weighted  to  more  closely  reflect  the  region  of  influence 
of  the  characteristic  lines  (refer  to  Figure  7) , eventually 
becoming  totally  hyperbolic. 


Admittedly,  such  a system  does  not  meet  all  the  criteria 
for  a potential  equation  solver  set  out  above,  but  represents  a 
workable  compromise  between  complete  agreement  with  the  above 
criterion,  and  the  practical  constraints  of  computer  code 
development . 


The  viscous  equation,  on  the  other  hand,  represents  a much 
easier  form  to  implement.  Dodge^  and  Dodge  and  Lieber2  have 
discussed  many  of  the  difficulties  associated  with  marchinq 
equation  solution.  For  this  program,  the  streamwise  diffusion 
term  was  determined  from  the  previous  iteration.  In  separated 
flow  regions,  the  difference  star  for  the  convection  term  was 
turned  upstream,  utilizing  the  values  established  on  a previous 
iteration.  With  these  adjustments,  the  solution  to  the  viscous 
equation  proceeds  as  a normal  marching  solution  to  a parabolic 
equation.  The  component  of  velocity  normal  to  the  wall  in  the 
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Figure  6.  Transonic  difference  stars. 
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region  between  the  wall  node  and  the  first  relaxation  node  was, 
however,  determined  from  the  continuity  equation.  This  was 
found  in  past  studies  by  Dodgel  to  produce  somewhat  more 
accurate  results. 

The  program  developed  to  implement  this  numerical  solution 
method  consisted  of  a modularized  code  written  in  a structured 
FORTRAN  language.  A basic  logic  diagram  for  the  program  may  be 
found  in  Figure  8. 
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Figure  8.  Program  logic  diagram 


SECTION  IV 
CALCULATION  RESULTS 


The  numerical  method  discussed  herein  has  been  utilized  to 
predict  separation  of  a two-dimensional  laminar  flat-plate  shock 
boundary  layer  interaction  flow.  The  test  case  used  was  that  of 
Hakkinen,  et  al7,  consisting  of  a laminar  flat-plate  flow  with 
an  upstream  Mach  number  TIq  = 2.0,  and  a Reynolds  number 
Re^  = 2.96  x 10^,  where 


W L 

o 


and 

W = freestream  inlet  velocity 
(WQ  = 1664  ft/sec) 

v = freestream  inlet  value  of  kine- 
matic viscosity  at  static  con- 
ditions 

L = distance  from  leading  edge  to 
shock  intersection  with  plate 
(L  = 0.1624  ft) 

An  oblique  shock  wave,  with  shock  angle  8 = 32.59°,  intersects 
the  flat  plate  as  indicated  in  Figure  9. 

The  grid  system  employed  in  the  numerical  solution  consisted 
of  a linear  relaxation  grid  and  a marching  grid  composed  of  the 
relaxation  nodes,  with  an  additional  nonlinear  grid  segment  near 
the  wall.  This  scheme  permits  a coarse  grid  to  be  used  for 
the  relaxation  solution,  and  yet  provides  a more  finely-spaced 
grid  within  the  boundary  layer  for  the  viscous  marching  equation 
solution.  Figure  10  illustrates  the  relative  spacing  of  the 
relaxation  and  marching  grid  nodes  within  the  boundary  layer. 

Results  of  the  present  numerical  method  were  compared  with 
both  the  experimental  data  of  Hakkinen,  and  results  of  the  time- 
dependent  scheme  of  MacCormack®. 

Figure  11  compares  surface  pressure  distribution  for  the 
three  sets  of  data.  Clearly,  the  pressure  plateau  in  the  region 
of  the  shock/plate  intersection  is  not  indicated  by  the  present 
method.  A much  sharper  pressure  rise  is  predicted  than 
evidenced  by  experiment. 
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As  a consequence  of  this  too-sudden  rise  in  static 
pressure,  the  predicted  values  of  skin  friction  coefficient, 

Cf , shown  in  Figure  12,  fail  to  decrease  as  far  upstream  as  the 
experimental  data.  The  predicted  region  of  reversed  flow  is 
also  much  smaller  than  that  indicated  by  Hakkinen's  data. 

The  predicted  values  of  Cf  are  presented  as  average  values 
within  a bandwidth  of  error  representing  the  standard  devia- 
tion of  Cf  as  calculated  at  as  many  as  8 nodes  within  the  con- 
stant pressure  and  density  portion  of  the  boundary  layer.  The 
skin  friction  coefficient  is  defined  by 
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P v 

W w 


1 
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p W ‘ 

OO  CO 
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where  the  subscript  w refers  to  condition  at  the  wall,  and  the 
subscript  00  refers  to  the  freestream.  A considerable  fluctua- 
tion in  values  of  3Wf/9xT  takes  place,  apparently  due  to  the 
nonlinear  character  of  tne  grid  system  near  the  plate. 

Figures  13  through  15  are  comparisons  of  boundary  layer 
velocity  profiles  between  the  present  steady-flow  method  and 
MacCormack's  technique.  In  general,  agreement  is  fair;  however, 
once  again  it  is  apparent  that  the  present  inethod  predicts  a 
smaller  region  of  separated  flow,  and  generally  larger  displace- 
ment thicknesses.  Discrepancies  in  the  assumed  inlet  boundary 
layer  profile  are  also  apparent.  Note  that  dimensions  in  the 
X2  direction  have  been  normalized  by  Hf  (where  Hf  = 0.003  ft)  — 
a grid  parameter  used  in  MacCormack's  data  presentation. 

Figure  16  presents  static  pressure  contours  over  the  solu- 
tion space,  based  upon  the  results  of  the  steady-flow  method 
discussed  here.  The  incident  and  reflected  shocks  are  apparent, 
along  with  a spurious  downstream  reflection  from  the  upper  bound- 
ary of  the  region.  The  effect  of  the  reflection  is  clearly 
visible  in  the  downstream  portion  of  the  skin  friction  and  sur- 
face pressure  curves. 
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Figure  13.  Velocity  profile  upstream  of  separation. 
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igure  16.  Static  pressure  contours 
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SECTION  V 

HOW  TO  IMPROVE  ACCURACY 

Because  of  time  and  fiscal  constraints,  the  development  of 
this  method  halted  with  the  results  shown  above.  It  is  the 
opinion  of  the  authors,  based  on  an  analysis  of  the  computed 
results  to  date,  that  considerable  improvements  in  accuracy  are 
possible  by  improving  the  accuracy  of  both  the  relaxation  and 
marching  numerical  methods.  The  relaxation  nodes  need  only  be  a 
subset  of  the  total  grid.  However,  relaxation  is  responsible  for 
satisfying  the  continuity  equation.  The  effects  of  the  increased 
displacement  thickness  downstream  of  the  shock  are  propagated 
forward  through  the  subsonic  region  of  the  flow  field,  creating 
the  shock  spreading  observed  above.  It  is  of  primary  importance 
that  this  effect  be  properly  accounted  for.  Figures  17  and 
13  show  the  results  of  the  above  calculation  on  the  mass  balance 
of  each  grid  element.  In  the  relaxation  grid,  agreement  is  very 
good;  but  inside  the  marching-only  region,  substantial  error 
occurs.  This  error  comes  about  from  several  sources  associated 
with  the  nonlinearity  between  the  last  relaxation  node  and  the 
wall.  A finite-difference  approach  to  numerical  method  deriva- 
tion assumes  that  a match  of  appropriate  derivatives  at  a point 
yields  a reasonable  estimate  of  the  founding  equation  in  the 
vicinity  of  the  match  point.  This  is  only  true  when  variations 
in  dependent  variables  are  nearly  linear,  or  at  most  quadratic. 

In  a boundary  layer,  such  as  shown  above,  this  assumption  is  in 
doubt.  Substantial  errors  can  occur  in  the  numerical  grid  trans- 
formation, the  evaluation  of  the  source  term  for  the  relaxation 
equation,  the  evaluation  of  derivatives  across  the  jump,  and  the 
like.  The  continuity  error  shown  in  Figures  17  and  18  is  the 
result  of  the  sum  total  of  these  effects.  Undoubtly,  a similar 
momentum  and  energy  error  also  occur.  The  cure  is  to  first 
produce  a marching  system  that  is  derived  from  an  integral 
equation.  Thus,  it  can  be  integrally  conservative  through  rapid 
changes  in  dependent  variables.  Secondly,  an  integral  average 
for  a source  term  of  the  potential  equation  must  be  established. 
As  an  example,  for  the  position  shown  in  Figure  17,  V • U,  part 
of  the  source  term,  differs  by  a factor  of  9 between  the  wall 
and  the  first  relaxation  node.  The  net  effects  of  the  source  are 
not  then  well  modeled  by  utilizing  its  value  at  the  relaxation 
node,  but  rather  a weighted  average  across  the  span  would  be 
needed.  A grid  system  that  is  more  nearly  uniform  would  also 
help.  It  is  obvious  when  examining  the  skin  friction  coeffi- 
cient, as  an  example,  that  for  W to  deviate  so  far  from  linear 
very  near  the  wall,  there  must  be  substantial  error  in  the  com- 
putation of  the  numerical  transformation  functions.  The  super 
position  of  two  fixed-spaced  grids,  as  utilized  by  MacCormack, 
would  substantially  alleviate  this  difficulty. 
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SECT I OH  VI 
CONCLUSIONS 

The  purpose  of  this  program  has  been  to  compare  a unique 
method  of  numerically  solving  the  steady-state  Navier-Stokes 
equations  with  experimental  data  unconfused  by  turbulent  model 
choice.  This  goal  has  been  met.  On  the  plus  side,  the  follow- 
ing can  be  stated: 

o The  method  has  met  its  rapid  speed 
potential,  with  run  times  in  the 
vicinity  of  1 minute  on  a CDC 
7600  for  a 3060  node  problem. 

Run  times  could  be  improved  sub- 
stantially by  more  careful  pro- 
gramming, and  better  selection  of 
the  grid  system. 

o A suitable  method  for  solving  the 

mixed  elliptic-hyperbolic  potential 
equation  has  been  developed. 

o Iteration  between  marching  and 
potential  equations  converged 
rapidly. 

o A simple  method  for  handling  both 

separation  and  diffusion  terms  within 
the  framework  of  a single-pass  marching 
equation  has  been  utilized  successfully. 

o Calculated  results  are  in  general  agree- 
ment with  experiment,  although  pressure 
rise  is  not  spread  over  as  great  a distance 
as  experimentally  observed. 

o A method  has  been  developed  that  could 
easily  be  extended  to  three-dimensions. 

Agreement  with  experimental  data  could,  however,  certainly 
stand  improvement.  It  is  expected  that  the  route  to  such  improve- 
ment is  through  a more  care':il  treatment  of  the  numerical 
details  involved  in  the  solution  technique,  and  does  not  reflect 
any  fundamental  limitation  of  the  method. 
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